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Abstract 

In Systems Biology there is a growing interest in the question, whether 
or not a given mathematical model can admit more than one steady state. 
As parameter values (like rate constants and total concentrations) are of- 
ten unknown or subject to a very high uncertainty due to measurement 
errors and and difficult experimental conditions, one is often interested in 
the question, whether or not a given mathematical model can, for some 
conceivable parameter vector, exhibit multistationarity at all. A partial 
answer to this question is given in Feinberg's deficiency one algorithm. 
This algorithm can decide about the existence of multistationarity by an- 
alyzing a, potentially large, set of systems of linear inequalities that are 
independent of parameter values. However, the deficiency one algorithm 
is limited to what its author calls regular deficiency one networics. 
Many realistic networks have a deficiency higher than one, thus the algo- 
rithm cannot be applied directly. In a previous publication it was sug- 
gested to analyze certain well defined subnetworks that are guaranteed 
to be of deficiency one. If these subnetworks are regular, then one can 
use the deficiency one algorithm to establish multistationarity. Realistic 
reaction networks, however, often lead to subnetworks that are irregular, 
especially if metabolic networks are considered. Here the special struc- 
ture of the subnetworks is used to derive conditions for multistationarity. 
These conditions are independent of the regularity conditions required by 
the deficiency one algorithm. Thus, in particular, these conditions are 
applicable to irregular subnetworks. 

1 Introduction 

In Systems Biology there is a growing interest in the question, whether or not 
a given mathematical model can admit more than one steady state. In cell 
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cycle regulation, for example, one can identify different phases of the cell cycle 
(Gl, S, G2 and M-phase) as different stable steady states. The cycle itself can 
then considered as a switching between these steady states. As parameter val- 
ues (like rate constants and total concentrations) are often unknown or subject 
to a very high uncertainty due to measurement errors and and difficult exper- 
imental conditions, one is often interested in the question, whether or not a 
given mathematical model can, for some conceivable parameter vector, exhibit 
multistationarity at all. 

A partial answer to this question is given in Feinberg's chemical reaction 
network theory, that links the ability of a mathematical model to exhibit mul- 
tistationarity to the structure of the underlying biochemical reaction network 
[3 [6l [3 [8] . The deficiency one algorithm, in particular, can decide about the 
existence of multistationarity by analyzing a, potentially large, set of systems 
of linear inequalities that depend on the network structure alone, that is, that 
are independent of parameter values. If any of these inequality systems is fea- 
sible, then multistationarity is guaranteed and one can compute steady states 
and rate constants from its solution set. If all are infeasible, then multista- 
tionarity is impossible, for any conceivable parameter vector (see, for example, 
[51 [5]). Observe that, in particular, this algorithm can also be used to prove 
that multistationarity is impossible. 

However, the deficiency one algorithm is limited to what its author calls 
regular deficiency one networks [5]. Many realistic networks have a defi- 
ciency higher than one, thus the algorithm cannot be applied directly. In [3l [9] 
we therefore suggested a way to circumvent this: instead of analyzing the com- 
plete network we propose to analyze certain well defined subnetworks that are 
guaranteed to be of deficiency one. If these subnetworks are regular, then one 
can use the deficiency one algorithm to establish multistationarity. If this is 
successful, then [3] gives sufficient conditions that are computationally simple 
to check to extend multistationarity from the subnetwork to the overall network. 

Realistic reaction networks, however, often lead to subnetworks that are 
irregular, especially if metabolic networks are considered (see e.g. [TU] for an 
analysis of the upper part of glycolisis). Consequently, the deficiency one algo- 
rithm cannot be applied to these subnetworks. If this irregularity is of a special 
kind (termed 0-irregularity in |10| ) , then one can regularize the subnetwork and 
apply the deficiency one algorithm to the resulting regularized subnetwork. It is 
then possible to extend multistationarity - so it exists - to the overall network 
using the aforementioned results of |3j. 

Here we follow a different approach: instead of trying to regularize a sub- 
network, wc use the special structure of the subnetworks defined in [3] to derive 
conditions for multistationarity. These conditions are independent of the regu- 
larity conditions required by the deficiency one algorithm. Thus, in particular, 
these conditions are applicable to irregular subnetworks. Of course it is still 
possible to use the results of [3] to extend multistationarity (once it can be 
established in the subnetwork). 
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2 Notation 



Consider the following (bio)chcmical reaction network with n 
B and with rn = 5 complexes A, 0, B, A + B and 2 A and r = 



: 2 species A and 
6 reactions: 



A 



ki 



k2 



k3 



k4 



k5 



2A 



Let X € M" be the vector of species concentrations (e.g. let xi be the concentra- 
tion of A and X2 be the concentration of B). By associating each concentration 
with the corresponding unit vector a of Euclidean space {A with ei and B with 
62 in case of the example) one can define m 'complex'-vectors yi (in case of the 
example yi = ei for A, j/2 = 0, the 2-dimensional zero vector for the complex 0, 
2/3 = 62 for B, 1/4 = ei + e-z for A + B and 2/5 = 2 ei for 2 A). Collect these in a 
matrix Y G 5?"^™. For the example one obtains 



y = [ 



yi ■■■ y5 \ 

10012 

00110 



Let la be the incidence matrix of the graph associated to the reaction network 
in standard form as defined in [7l [8], that is a graph, where node labels are 



unique. This means that one has la £ { — 1,0 
obtains 

^ -1 



la = 



ly 





-1 
1 



For the example one 



Finally let k G IR^q be the vector of rate constants, that is for the example: 

fc = (fci, . . . , fcg) ■ 

The stoichiometric matrix N is defined as the product 

N:^YIa (1) 

of the matrix of complexes Y and the incidence matrix of the associated directed 
graph la. 

Definition 1 (Reactant Complex, Educt, rh). A complex that has at least one 
outgoing edge is called a reactant complex. We use the symbol fh < m to denote 
the number of reactant complexes. 

Let y be a reactant complex. Then all species with indices contained in supp (y) 
are called educts. 
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For simplicity we assume - w.l.o.g. - the following ordering of complexes: 

Assumption 2 (Complex Ordering). Assume that the complexes are ordered 
such that the first rh complexes are reactant complexes. 

Under this assumption the mapping reac that associates every reaction with 
its reactant complex has a particular simple form: 

Definition 3 (Mapping reac). Let reac : {1, r} {1, m} be defined 
as 

reac (j) = i, yi is the tail of reaction j. (2) 

If mass-action kinetics is used, then the reaction rate Vi {k,x) associated to 
the i-th reaction is given as the monomial Vi (fc, x) = ki x^""-"'"' (i.e. the reaction 
rate Vi {k,x) is proportional to the product of the educt concentrations). One 
obtains the following function v(k, x): 

Definition 4 {v{k, x), $(a;), (x)). Using mass action kinetics, the vector of 
reaction rates is defined as 

v{k, x) := diag (fc) $ (x) , (3a) 

where 

$ (x) := (x^— , . . . , x?'— (-) f . (3b) 
Let Ci denote the unit vectors of Euclidian n-space and define 



n := 



T 

^rcac(l) 



T 

rcac(r) 



(3c) 



*(^) :=(^^'W. (3d) 
Note that this implies that 

$ (.t) = n ^' (x) (3e) 

and thus 

v{k, x) = diag (fc) $ (x) = diag (fc) n 5' {x) (3f) 

hold. Let 

^ = [y.W...,^ (3g) 

he a matrix having the exponents of"^ [x) as column vectors (recall assumption\^ 
and note that this implies that Y contains the first ffi columns ofY). 
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For the example one obtains 



v{k, x) = (fcixi, k2, ^3, k4X2, k5,xiX2, kexl) . 
Observe that Y ~ Y and thus 

$ (x) = * (x) , 

in this case. Then the following system of Ordinary Differential Equations 
(ODEs) describes the dynamics of the species concentrations: 

X = Nv{k, x), (4a) 

If the stoichiometric matrix N E J?"^'' does not have full row rank, the system 
is subject to 'conservation relations': let s = rank(iV) < n, then there is a 
matrix W G jR"x»-'' with W'^ N = and 

x{t) = c (4b) 

along solutions x{t) of l|ia)) . cf. [T]. 

As we are mainly interested in positive steady states, the pointed polyhedral 
cone ker(y /q) n M^q is of particular interest. The symbol E is used to denote 
the unique (up to scalar multiphcation) generators of ker(y la) Ci IR^q. Let p be 
the number of generators, if p > 1, then E G FT^^ is a matrix whose columns 
are the generators of ker(y /„) n J??>q, if p = 1, then E G ST is a vector. Define 
the set of all nonnegative vectors x G J^>o such that E x \s positive: 

A(£:) := {x G l?|o I > o} . (5) 



3 Some remarks about positive steady states 

The structure of (|4a|) motivates the following result concerning positive steady 
states: 

Lemma 1 (Existence of positive steady states). Consider a system of ODEs 
as in |.^ap , with stoichiometric matrix N and let E G Sl^^^ be the generator 
matrix o/ker(y /a)n-/R>g. Let k G H^q be given. Then the positive vector a is 
a solution to the polynomial equation N v{k, a) = 0, if and only if there exists 
a vector A G A {E) with 

/s = diag($(a-i)) E\. (6) 

Proof. Follows from the fact that a > and A: > implies f (fc, a) > 0. Thus 
Nv{k, a) = holds if and only if v{k, a) G ker(y /„) DM^q, that is, if and only 
if v{k, a) = EX, for some A G A(i?). As v{k, a) = diag (A:) cf> (a) ^ follows 
immediately. □ 
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Remark 1. If a positive steady state exists, then must hold. The condition 
0j can thus be used to constrain the set of rate constants that allow the existence 
(of at least one) positive steady state. 

Remark 2 (Positive steady states) . Consider a system of ODEs as in |^a[ ) and 
let E G JR^yi)' be the generator matrix o/ker(y/a) H Hy^- The system has a 
positive steady state, iff her {Y la) H ^ and the rows of E are nonzero. 

Remark 3. Consider a system of ODEs as in J-^ap and let E G -2?>g^ be the 
generator matrix ofkei^VIa) H IR^q and suppose that E does not contain any 
zero rows. Then every positive vector a can be a steady state of |4a[ ), by choosing 
k as in where A G A (E) is free and takes the role of the rate constants. 

4 Subnetworks defined by stoichiometric gener- 
ators 

In this section the following concepts from graph theory will be used (two of 
them are standard definitions in graph theory, that stated here merely for con- 
venience, while the third, very common in CRNT, is derived from those two): 

Definition 5. !lllf Connected component: the maximal connected subgraphs 
of a graph 

Strongly connected component: a directed graph is called strongly connected 
if there is a path from each vertex in the graph to every other vertex. The 
strongly connected components (SCC) of a directed graph are its maximal 
strongly connected subgraphs. 

JBjl Terminal strongly connected component: an sec that has no outgoing edge 

Next we recall some results concerning subnetworks defined by stoichiometric 
generators 

Lemma 2 (Properties of subnetworks defined by stoichiometric generators). 

For a subnetwork that is defined by a stoichiometric generator E the following 
properties hold: 

(a) Graph of the network in normal form is a forest of trees 

(b) Terminal strongly connected components consist of a single node (complex) 

(c) The deficiency of the network is one 

(d) The deficiency of every connected component is zero 

(e) If every connected component contains only one terminal strongly con- 
nected component, then the network is regular (in the sense of CRNT, 

example) 
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(f) kei{N) = [E] 
Proof. 

(a),(b) Follow from the definition of the generators of ker {Y la) H 
(c)-(f) A proof can be found in [3] 

□ 

The following corollary is an immediate consequence of Lemma [21 Q and 
Remark [3] 

Corollary 1. Consider a biochemical reaction network that is defined by a 
stoichiometric generator. Then 

(i) any positive vector a is a steady state of \4a\ l, if k is chosen as in @) 

(ii) for an arbitrary but fixed positive a, k as in (0j is fixed up to scalar mul- 
tiplication (i.e. the positive X) 

(Hi) (positive) scalar multiplication of k corresponds to a time scaling of the 
ODEs, thus one can - w.l.o.g. - choose A = 1 

From here on we assume that the system has at least one positive steady 
state, that is 

Assumption 6. The vector of rate constants is given by 

fc = diag($(a-i)) £; (7) 

for some a g J^>o- 

Consider Lemma [T] and especially the facts that for networks defined by 
stoichiometric generators E consists of one (column) vector and that - w.l.o.g. 
- A = 1. Then the ODEs (|4ap arc equivalent to 

X ^ Nv{k, x)=:N diag (E) diag ($ {a-^)) $ (x) = N diag (E) $ (^^) , 

where $ (a:) = (x) (cf. Definition [4]) . Thus 

x^N diag (E) n * = iV diag (E) U diag (* (a"^)) * (x) (8) 

follows. 

Remark 4. Systems like ^ are sometimes called generalized mass action sys- 
tems. For those systems reaction rates Vi {k,x) are still defined as monomials 
ki a;^"'"=(') , however the exponent vector yrcac(i) does not need to correspond to 
the reactant stoichiometry anymore. 

Further observe that for the special system defined in (0j '5 (a^^) takes the 
role of the rate constants. 
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To establish multistationarity we need to show the existence of a second 
steady state b G JR"q with 

Nv{k, b) = 0, 

for the same vector k. That is b must satisfy: 

N diag {E) n diag (* (fl-^)) (6) = (9) 

Obviously ker [N diag [E) 11) = [l_p] (as N is the stoichiometric matrix of a 
subnetwork defined by a stoichiometric generator; to sec this recall that (i) 
ker [N) ~ [E] and (ii) 11 has full column rank and (iii) row vectors of 11 are unit 
vectors: thus diag (i?) 111^ ^ E). It follows that Q is equivalent to (observe 
that a, 6 > implies * (^) > 0): 




= aim," > 



Apply In (•) to obtain the linear system 

y^M = In («)!,-„, (10a) 

where ^ 

;,:=ln- = fln^,...,ln^^ . (iQb) 
a \ ai an J 

The previous discussion motivates the following Lemma: 

Lemma 3 (Parameterizing positive steady state solutions). Consider the ODEs 
derived from a biochemical reaction network that is defined by a stoichiometric 
generator. Let 

:= |m e 5?" I 3p > 0, such that f'^^ = pl„-j| . (11a) 

IJ M- ^ 0, then fi E A4 and a G Ji!>o parameterize positive solutions of the 
polynomial equation N v{k, b) = 0. Let jjL E M and let 

fc = A$ (a"i) , A > (lib) 

be the vector of rate constants. Further let 

6 = diag (e^) a. (11c) 
Then Nv{k, a) and N v{k, b) = hold. 

Proof. Let /i G let k and b be as in (jllbp . picp . respectively. Observe that 
Nv{k, a) = follows from Lemma[TJ Thus we have to show that Nv{k, b) — 
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holds. To this end observe that 

Nv{k, b) ^ N diag{XE) diag ($ (a"^)) $ (&) 

= XN diag (E) diag ($ (a^^)) $ (diag (e'') a) 
^ \ N diag (i;) diag ($ (a^^)) diag ($ (a)) $ (e'^) 
= AA^ diag(i;) n^'(e^) 

^ \ N diag (S) n e"^"^ 

= A7V diag(i;) ne''l„^ = 0. 



□ 



Remark 5. If either (i) If^ G 



or ('iij ker {^"^^ is nontrivial (i.e. rank < 
n), or both, then a fixed vector a € IR^o together with k as in Hi 1 b\) defines an 
infinite set of positive steady states b G J^>o- To see this assume that (i), (ii) 
or both hold. Then \10a\) is solvable and the solution set 



M := |/i e I 3/9 > 0, such that Y"^ = p 1„-,| 

defines a linear subspace, that is there exists a matrix M and a vector k of 
appropriate dimensions, such that 

fi e M ^ f-i = M K 

holds. As every fi £ M defines a positive b and as a linear vector space contains 
infinitely many elements fi, there exists, for a fixed a G ^?>o infinitely many 
b e iR"o with 



ba (k) = diag(e^^'=) 



5 Subnetwork multistationarity 

If the conditions of Remark O hold, then there exists an infinite set of posi- 
tive steady states, even for a given a G -ff?"g (recall that 5* (a~^) takes the 
role of the rate constants). Fix a G iR>o and thus k = diag ($ (a^^)) E = 
diag (i?) n^' (a^^)- Then all ba (k) as defined in are steady states. How- 
ever, for a given initial condition xq G -K"q the system 'sees' only a subset of 
set of positive steady states. 

To see this recall the ODEs derived from a biochemical reaction network 
X = Nv{k, x) and let s := rank (A^) < n. Let S, W be orthonormal bases of 
[N] =: S and S'^, respectively. Similar to [T], introduce the transformation 

y — S'^ X, z — X and x = £^{x,y) = S y + W z. 

In the new coordinates the ODEs read as 

y = 5^ i = 5^ iV diag [E) n diag (a-^)) * {y, z)) 
z = w'^x = 
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showing the hivariance of S. Let x (0) = xo £ K>oi then the sohition x {t) is 
given by 

x{t) = xo+ [ N diag (E) U diag (* (a~^)) * (a; (r)) dr. 

"'0 

For the new coordinates one obtains 

y (t) = S^xo + S^ f N diag {E) U diag (a"^)) * (x (t)) dr 
Jo 

z (t) = W'^ xq = const. 

that is solutions are confined to paraUel translates of S. Thus, for a given initial 
condition, the system 'sees' only those positive steady states that are in the 
intersection of ba (k) and S. Observe that S := [N diag (E) 11] C 5 and that, 
by a similar argument, solutions are confined to parallel translates of S. This 
motivates the following definition of multistationarity with respect to a linear 
subspace as introduced in [§]: 

Definition 7. Given a subspace V C IR", the system x — Nv{k,x) from ^a\ l 
with stoichiometric subspace 5 = im {N) is said to exhibit V -multistationarity 
if and only if there exist a positive vector k 
positive vectors a, b £ -K"q with 



G J^>o ^''T-d o,t least two distinct 



Nv{k,a) =0, 



Nv{k,b) =0, 
b-aeV. 



(13a) 
(13b) 



Remark 6. Note that ifV = S, then Definition^ is equivalent to the familiar 
definition of multistationarity in Chemical Engineering and especially in CRNT 
as defined, for example, in fH 

Remark 7. Note that for subnetworks defined by stoichiometric generators two 
linear subspaces are of particular interest: [N diag (E) H] and the image of 



stoichiometric matrix of the overall network 



N 



Multistationarity with respect 



to [N diag (E) 11] means that the subnetwork can exhibit multistationarity, if it 
is considered in isolation, while multistationarity with respect to N 
that the subnetwork as part of the larger network can exhibit multistationarity. 
Thus, if a subnetwork exhibits N -multistationarity, but not [N diag (i?) 11]- 
multistationarity, then this subnetwork can give rise to multistationarity for the 
overall network, even though, in isolation, it does not exhibit multistationarity. 



Remark 8. As an illustration consider the network in Fig 1(a) For this net- 
work the matrix N is given by 



N = 



10-1 
1-1 
1-1 



(14a) 
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the unique generator of her {Y la) H -/R>o is given by 

E={1, 1, 1, if. 
For $ (x) one obtains 

= (1, 1, X3, XiX2f 

and therefore 



(14b) 



(14c) 



n = 



1 

1 

1 

1 



^ [x) = (1, X3, Xi X2f 



(14d) 



Thus one obtains 



1, 1, 



ai a2 as 



(Me) 



for arbitrary a > 0. Observe that for this example N diag (E) 11 = TV 11. It is 
straightforward to verify that all points on the following one- dimensional curve 
are steady states (parameterized by p > 0): 



1 

P, 1 
P 



(14f) 



As N has full row rank, the left kernel is trivial, that is, it is spanned by W = 0, 
the three-dimensional zero vector. Pick two distinct real numbers pi >, p2 > 0- 
Then Xg [pi), Xg {P2) O'fc steady states that satisfy {xs (pi) — Xs {P2)) — 
(i.e. . Xs (pi) — Xs {P2) G [N]. Thus, according to our definition, the system 
exhibits [N]-multistationarity. 

It fails to exhibit [N H] -multistationarity, as for a given initial condition 
Xq > 0, all trajectories converge to a unique steady state. That is due to the 
fact that 

"1-1 

1 -1 
1 -1 



NU = 



has a nontrivial left kernel W = (1, — 1, 0) . Thus, for a given initial condi- 
tion, a trajectory is confined to an affine linear subspace perpendicular to W . 
And all trajectories starting in particular affine linear subspace converge to the 
same steady state (demonstrated numerically for a = 1 in Fig 1(c) and 1(d)). 
Note that, using k as in {l-jc^ the ODEs are equivalent to a system of ODEs 



derived from the network displayed in Fig 1(b) . a weakly reversible deficiency 
zero network. From the Deficiency Zero Theorem ^ follows that this network 
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k4 kg 

A + B ^0 



(a) Example network 




B 



A + B 



(b) Transformed network 




T 1 1 r- 




7 8 9 10 



(c) Simulation for selected initial conditions 



(d) Projection in the x'i-2;2-plane 



has a unique, asymptotically stable positive steady state - relative to a given 
initial condition. 



Further note that the network in Fig. 1(a) is irregular in the sense of CRNT 



(cf. l^) and that it fails to exhibit [N Ii\-multistationarity, while it exhibits [N]- 
multistationarity. 



6 Establishing multistationarity for subnetworks 

Consider a biochemical reaction network defined by a stoichiometric generator. 
From Section m it is known, that for a given a S -K>o all points ba (k) as defined 
in (fT2|) are steady states. Moreover, the set a linear subspace, as defined 
in (jllap contains all yu = \nba{n) — In a. From Section [5] it is known, that 
V-multistationarity requires 



In 6a (k) -\na <E M 

To this end a result from [2] can be used. To state the result, let sign [u] denote 
the sign pattern of the vector u € . Then v = sign (u) is a vector with entries 
Vi € {+, — , 0} depending on whether Ui > 0, m < or Ui = 0, respectively. 
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Lemma 4 (cf. [5]). Let Mi C JR" and M2 C JR" &e iwo nontrivial subsets of 
IR"- and define M3 := { (mi, TO2) £ -^1 x M2I sign (mi) = sign (m2) } as the set 
of all ordered pairs (mi, m2) of elements mi G Mi and m2 £ M2 wii/i ifte same 
sign pattern. Two positive vectors p and q with hiq — Inp G Mi and q — p E AI2 
exist, if and only if M3 7^ 0. Then p and q are given by 



\Pi> 0, if mii = 0, 
where pi denotes an arbitrary positive number and 

(90.=i,...,„=e™-p.. (16) 

Thus - using ba (k) instead of q, a instead of p and Ai instead of Mi, V 
instead of M2 - all one has to do is to find a vector ji € M. and a vector s G V 
with sign (/i) = sign (s) to establish V-multistationarity. 
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